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ABSTRACT 

The linearization of the meteorological equations around a specified reference state, 
usually applied in NWP to define the linear system of constant-coefficients semi-implicit 
schemes, is outlined as an unnecessarily restrictive approach which may be detrimental in 
terms of stability. It is shown theoretically that an increased robustness can sometimes be 
obtained by choosing the reference linear system in a wider set of possibilities. The poten- 
tial benefits of this new approach are illustrated in two simple examples. The advantage 
in robustness is not obtained at the price of an increased error or complexity. 
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1 Introduction 



The semi-implicit (SI) technique was proposed in the 70's (Robert et al, 1972) as a 
suitable and efficient method for solving numerically the partial differential equations used 
in meteorology. At this time, the SI technique was applied to hydrostatic primitive equa- 
tions (HPE), and its success in this context made it very popular in the field of numerical 
weather prediction (NWP). The suitability of the SI technique for the fully elastic Euler 
equations (EE) was then advocated (Tanguay et al. 1990), with the aim of combining the 
advantages of a system valid at any scale and an efficient time-discretization, as required 
for NWP purposes. 

The essence of SI schemes is a linear separation of the source terms of the complete 
system to be solved, with an implicit treatment of this linear part. For the purpose of 
this paper, three main types of SI schemes can be distinguished. The coefficients of the 
implictly-treated linear terms can be : (i) constant in time and horizontally homogeneous 
(Simmons and Temperton, 1997; Bubnova et al., 1995, Caya and Laprise, 1999); (ii) 
constant in time only (Thomas et al., 1998; Qian et al., 1998); and (iii) non-constant 
(Skamarock et al., 1997, Cullen et al., 1997). 

This paper only considers SI schemes belonging to the class (i), which are designed 
under "constant-coefficient SI schemes" in the following. However, it should be outlined 
that since only the separation of thermal terms is considered, all results and conclusions 
extend identically to those SI schemes of type (ii) for which the reference temperature is 
horizontally homogeneous (e.g. Thomas et al., 1998; Qian et al., 1998). 

The underlying principles usually applied in the design of constant-coefficients SI 
schemes are the following : 

(i) define a stationary SI reference basic state X* ; 
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(ii) linearize the meteorological system M. to be solved around this steady state, to 
obtain a linear system C* ; 

(iii) treat the linear part of the evolution C* with a centred-implicit scheme, and the 
remaining "non-linear" part (M. — £*) with a centred-explicit scheme. 

However, due to the explicit treatment of the non-linear (NL) residuals, the stability of 
this type of scheme is not formally guaranteed, especially with long time-steps. Indeed, the 
application of the above technique sometimes leads to unexpected unstable behaviours. 
The two following problems (referred to as PI and P2 hereafter) illustrate the kind of 
limitations which can be encountered with constant-coefficients SI schemes designed using 
the principles (i)-(iii) : 

PI : With HPE, the introduction of a vertically varying reference thermal profile T* close 
to the actual atmospheric profile, although reducing the magnitude of the thermal 
NL residuals, leads to a scheme which is less robust than when a warm isothermal 
T* profile is used (see e.g. Simmons et al. 1978, SHB78 hereafter). 

P2 : For two time-levels (2-TL) SI discretizations, the EE system is extremely unstable 
while the HPE system is stable, as discussed in Benard (2003, B03 hereafter). 

As mentioned above, the constant-coefficients SI technique has traditionally been ap- 
plied to NWP by explicitly following the three principles (i)-(iii), but this method is 
unnecessarily restrictive. 

As stated in B03, the SI scheme can be viewed as the very first iteration of a generalized 

pre-conditioned fixed-point algorithm for iteratively approaching the pure centred-implicit 

scheme. In this light, C* appears to be nothing else than the linear pre-conditioner of the 

fixed-point algorithm (this pre-conditioner is necessary in such an algorithm for allowing 

the convergence of the iterative process) . This point of view outlines the arbitrariness of 
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the choice of the C* system, provided a satisfactory convergence for the iterative algorithm 
is ensured. 

When facing unexpected problems as (P1)-(P2), a possible solution, advocated in 
this paper, is to relax the constraints (i)-(ii) and to seek £* deliberately as an arbitrary 
constant-coefficients linear system, i.e. not obtained through the linearization of M. around 
any reference state. This method is illustrated in the two following practical examples. 

2 Proposed solution to the problem (PI) 

It is a well-documented fact that if the stability of the SI scheme is obtained by forcing 

NL residuals to large values (in such a way that their sign is controlled), then the response 

of the scheme is deteriorated, especially from the point of view of phase-speed errors. 

When exaggerated, this strategy has a negative impact even on slower transient processes, 

making it unattractive for NWP. A natural way to alleviate this risk with certainty is thus 

to reduce the magnitude of NL residuals. This is precisely the idea which was tested in 

SHB78, by comparing the properties of SI schemes obtained when choosing isothermal 

and non-isothermal profiles of the reference temperature T*. The non-isothermal profiles 

were chosen close to standard atmospheric profiles, in such a way that the magnitude of 

thermal NL residuals was reduced compared to the case with isothermal T*. However, 

the experimental results in terms of stability were clearly worse for the non-isothermal 

option : when the tropopause of the T* profile was above the actual one, the scheme 

became highly unstable. Given this experimental fact, the recommended solution, widely 

followed afterwards, was to use a warm isothermal profile for T*, thus implicitly accepting 

to sacrifice a better response of the scheme for an increased robustness. 
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2.1 Analysis of SHB78 situation 

In this section, the HPE system in o coordinate is considered with a three time-level 
(3-TL) leap-frog SI time-discretisation. The theoretical framework proposed in B03 is 
used here to perform a stability analysis, and the reader is referred to this paper for more 
details on notations and algebraic developments. The framework is idealized in order to 
allow simpler analyses (Cartesian vertical (x, z) plane without orography ; dry, adiabatic, 
non-rotating equations). A resting state X with a thermal profile T(a) is considered. All 
atmospheric evolutions are assumed to consist in small perturbations around X (referred 
to as the "actual" state hereafter), and the meteorological system M. is linearized around 
this state, in order to allow tractable analyses. In the notations of B03, the system thus 
writes : 



?R -RQ— - RT d2q (1) 
dt dx 2 dx 2 



dT RT 



dt C p 



(2) 



where D is the horizontal wind divergence, T the perturbation temperature, q = ln(7r s ) 
, and ii s is the surface pressure. Note that in (J2J), the last RHS term is the contribution 
of vertical advections for T. This Eulerian form can be shown, in this linear framework, 
to be also valid for a semi-Lagrangian discretization, under the assumptions of perfect 
solution for the displacement equation and perfect interpolators, consistently with the 
current space-continuous context. A modified version of the cr-coordinate static-stability 
for the actual state X is introduced through : 



RT dT 

^=!T v - a Ta (4) 
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For the reference state X* used to define the SI scheme, a profile T*(cr) is also assumed, 
and the system is linearized around this reference state, according to the above principles 
(i)-(ii). The C* system and the static-stability 7* obtained through this procedure are thus 
formally identical to (0)-(JBl) and (pQ), respectively, simply substituting T* for T everywhere. 

These two static-stabilities (7, 7*) are now assumed uniform in the whole domain for 
the purposes of the analysis. A "non-linearity" factor is defined by a = (7 — 7*)/7*. It 
should be noted that the case of an isothermal SI reference state is also included in this 
formalism since it results in a uniform static-stability 7*. Following exactly the method 
presented in B03, the system (jH^fjBJ) is first transformed into an unbounded system : 



d \ dD 



a- 



da dt 



d\dT 



RV 2 T 



(5) 
(6) 
(7) 



The normal modes of the system are then : 



i/j(x,a) = $ exp(ikx)a {iu ~ 1/2) (8) 

(9) 

where (k, v) e R and ip represents either D or T. Pursuing the analysis as in B03, it is 
finally found that in the limit of long time-steps, the 3-TL SI scheme is stable for : 

0<7<2 7 *. (10) 

This result extends (and is fully consistent with) those obtained in previous related studies 

(SHB78, and Cote et al. 1983, CBS83 hereafter). Moreover, it allows an understanding 
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of the instability observed in SHB78 for their SI scheme with non-isothermal reference 
profiles : when the tropopause of the SI reference state is higher than the tropopause of 
the actual state, the above criterion ([TO]) is locally violated between the two tropopauses, 
resulting in an unstable scheme. However, for warm isothermal profiles of T*, the latter 
instability disappears, as empirically found by SHB78, because in this case, 7* has a 
high value at any level, and therefore is larger than 7/2 in all the depth of atmosphere, 
whatever may be the location of the actual tropopause. 

2.2 Proposed modification 

The fundamental difference between the two options examined by SHB78 is not in 
the values of T* themselves (which actually deviate marginally between the two types 
of considered reference thermal profiles T*), but in the presence or not of the advective 
term (dT* /da) in the C* system, because this term dramatically modifies the apparent 
static-stability 7*, as seen in (J1J. 

Hence, according to the new approach proposed in this paper, a natural solution to 
ensure a more stable scheme while keeping a non-isothermal T* profile, is to deliberately 
remove the resulting advective term in the initial C* system. This modification can be 
expected to combine both the advantages of small residuals (because T* can be made 
closer to actual atmospheric thermal profiles) and optimum stability (because the ap- 
parent static-stability in the C* system is large at any level). It is worth noting that the 
mathematical structure of the C* system in this modified SI scheme with a non-isothermal 
T* profile is exactly the same as for a traditional SI scheme with an isothermal T*, hence 
the modification in any pre-existing application is straightforward. 

In order to illustrate the consequences of this modification, a situation close to the one 

examined in SHB78 is considered. A class of vertical thermal profiles is introduced by : 
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where T = 220 K, 70 = 30 K, and or is a varying parameter specifying the level of the 
tropopause. The value <x£ = 0.25 is chosen for the SI reference state T*, while for the 
actual state T, the tropopause level Wr is left as a free parameter in the interval [0.1, 0.5]. 
The static-stability is 7 = (R/C P )T = 62.9 K in the isothermal (T = T ) "stratosphere", 
and 70 = 30 K in the "troposphere" for both T and T* profiles. The only difference 
between the two piecewise-constant profiles of static-stability (7, 7*) is thus the location 
of their tropopause. 

The stability analysis is not straightforward for such multi-layers systems, hence the 
stability of the systems is diagnosed through vertically-discretized analyses exactly as 
in CBS83. In this method, the whole vertically- and time-discretized system for a given 
horizontal mode is considered as a linear "amplification matrix" acting on a generalized 
vertical state- vector, and the growth rate T of the system is the maximum modulus of the 
set of eigenvalues of the amplification matrix. The vertical structure of the most-unstable 
mode is given by the associated complex eigenvector. The vertical discretisation in the 
analyses presented here is the same as in Simmons and Burridge (1981), and is equivalent 
to the one used in SHB78. The vertical domain is described through 80 regularly-spaced a 
levels, and the analyses are performed for a mode with k = 0.0005m -1 (the results are not 
qualitatively sensitive to k). As dicussed in B03, the examination of the stability in the 
limit of long time-steps is relevant since long time-steps area target in NWP. The value 
chosen here is At = 2000 s (here also, smaller time-steps do not change qualitatively the 
conclusions). 

The growth-rates for the traditional SI scheme and for the proposed modified SI scheme 

are depicted as a function of Wr in Fig. [TJ For the traditional SI scheme, the results are 
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(11) 



fully consistent with the criterion obtained through the above analysis (it has also been 
checked that a slight increase of the tropospheric static-stability 70 from 30 to 35 K 
results in a stable scheme for any value of or in the explored interval, in agreement 
with the stability criterion derived above). For the traditional design used in SHB78, the 
SI scheme is unstable as soon as the actual tropopause is lower than its SI reference 
counterpart. For the modified SI scheme proposed here, the stability is obtained even in 
the previously unstable situation, and is thus comparable to the case with an isothermal 
T* profile. 

In this analytical context, the modified scheme reaches the initial aim of reducing 
the magnitude of NL residuals while ensuring a robust scheme. The extension of these 
theoretical results to fully realistic frameworks has not been investigated further, but the 
approach seems worth considering since it potentially combines the two advantages of 
robustness and accuracy. 

3 Proposed solution to the problem (P2) 

A stability analysis of the EE system in the space-continuous SI framework, proposed 
in B03, shows that the 2-TL time-discretisation is very unstable in the presence of thermal 
NL residuals, while the HPE system are acceptably stable in the same context. 

As in the above section, a theoretical analysis reveals the causes of this dramatic 
destabilization in simplified contexts. In the following, it is shown that the destabilization 
originates from the fact that the thermal NL residuals corresponding to the terms respon- 
sible for gravity and elastic waves systematically have opposite signs. To better illustrate 
this explanation, the following couple of excerpts from the complete linearized EE system 

in <j vertical coordinate [see (52)-(56) in B03] can be examined : 
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dD 
~dt 
dT 
~dt 



RQV 2 T 



RT 



D 



(12) 



(13) 



and : 



dd 

dt 
dV 
~dt 




(14) 



(15) 



where all notations follow B03. 

The first sub-system describes the horizontal propagation of gravity waves, while the 
second describes the vertical propagation of elastic waves. Neglecting the Boussinesq effect 
represented by the term "+1" in the RHS of JSJ), these two sub-systems are formally 
identical, the only noticeable difference being the location of the T factors (at numerator 
vs. denominator). As a consequence, for a given set of actual and reference temperatures 
(T, T*) the thermal NL residuals always have an opposite sign in the two systems. For the 
purposes of the analysis, the thermal profiles T and T* can be considered as isothermal. 
Let a = (T — T*)/T* be the thermal non-linearity parameter for the considered simplified 
problem. The stability properties of the first sub-system for a are thus the same as those 
of the second sub-system for — a/(l + a). Since the stability of the first sub-system (in 
2-TL SI) for long time-steps implies a < (see B03 for details), the stability of the second 
sub-system necessarily implies a > 0, and the stability domain for a complete SI system 
which would include the two types of waves thus vanishes. In other words, if T* is chosen so 
as to stabilize horizontally propagating gravity waves, then vertically propagating elastic 
waves will be unstable, and vice versa. The problem is of course not present for HPE since 
this system does not allow the propagation of elastic waves. 
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A natural solution to restore systematically the same sign for thermal NL residuals 
in the above two sub-systems, is thus to introduce different values of T* for each sub- 
system, that is : T* for the gravity- wave system (fT2|) - (fT3|) . and TJ. for the elastic- wave 
system fTH)- (fT5|) . Noting T% = rT*, the stability domains for the first and second system 
become a < and (r—1) < a respectively. As a consequence, choosing r < 1 allows a non- 
empty stability domain for long time-steps to be restored. In terms of temperature, the 
stability is then ensured if T B < T < T* in this isothermal context. The stability domain 
for T can thus be arbitrarily extended, by setting T* arbitrarily warm, and arbitrarily 
cold, with the limitation that an exaggeration in this direction finally deteriorates the 
response of the scheme, as outlined above. 

The application of this solution to the complete EE system is straightforward : for all 
occurences of T* at numerator in the initial linear system [i.e. the reference system £.* 
associated to (52)-(56) in B03], the traditional warm value T* should be kept, while for 
the occurences of T* at the denominator, the cold value Tg should be imposed. Here also, 
the modification from any pre-existing application is straightforward. 

The theoretical impact of this modification is first illustrated with a stability analysis 
of the complete EE system for 2-TL SI schemes in the context of isothermal T and T* 
profiles and linear evolution around X as in B03. The analysis for T E = T* is given in B03, 
and can be repeated in a formally similar way for the modified SI scheme T% ^ T* . The 
growth rates obtained in the long time-step limit for the initial and modified SI schemes 
are depicted in Fig. El for r = 1 and r = 0.5 (i.e. T£ = T*/2). The results are fully 
consistent with the above simple analyses : the modified SI scheme is found to be stable 
in the interval (r — 1) < a < 0, while the traditional SI scheme is always unstable. 

In order to evaluate the potential benefit of the proposed approach for NWP, the modi- 
fication was then tested in real-case conditions with the adiabatic semi-Lagrangian version 
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of the Aladin-NH model (Bubnova et al. 1995), used with a 2-TL SI time-discretization. 
The model was integrated for 3 hours for a randomly-chosen situation consisting of a 
strong flow over real topography, in a domain which includes the montanous Pyrenees 
region. The horizonal resolution is 2.5 km in horizontal directions, and the time-step is 
80 s. The vertical coordinate is the mass-based hybrid coordinate defined in Simmons and 
Burridge (1981), and the domain is discretised along 41 irregular layers with a thickness 
increasing with height, in the usual NWP fashion. Integrations are performed without any 
time-filter (see B03 for a discussion on the detrimental effects of time-filters in 2-TL SI 
EE system) . A weak fourth-order horizontal diffusion is applied to avoid the accumulation 
of energy in the smallest resolved scales during the course of the integration. 

Fig. [HI shows the evolution of the whole domain spectral norms of the horizontal vor- 
ticity ( and divergence D for the traditional and modified versions of the 2-TL SI scheme. 
The traditional SI scheme is used with T* = T|,=300K, and the modified SI scheme with 
T*=300 K, Tg=150 K. The original 2-TL SI scheme is clearly unstable, since the integra- 
tions diverge after 11 time-steps, while the modified 2-TL SI scheme behaves stably during 
the 3 hours of the integration. This experiment clearly indicates a potential advantage of 
using the modified SI scheme in NWP with 2-TL EE systems. 

4 Comments and conclusion 

All the discussions in this paper apply equally to 2-TL and 3-TL schemes. They can 
also be extended straightforwardly from SI schemes to the emerging class of iterated 
centred-implicit (ICI) schemes, as examined in B03, because these schemes are based on 
the same kind of linear separation of the meteorological system to be solved implicitly. 

For the problem (P2), the proposed solution offers a smaller interest for 3-TL schemes 
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than for 2-TL schemes, because 3-TL schemes already have a degree of robustness com- 
patible with a NWP use as far as thermal NL residuals are concerned. Nevertheless, the 
proposed solution is believed to be worth considering for high resolution modelling with 
the EE system in combination with 2-TL schemes For the EE system discretized with 
2-TL SI schemes, it may allow in particular to remove the strong time-filters used so far, 
with their detrimental effects in terms of response. Moreover, it would be interesting to 
extend this modification to systems in height-based coordinates. It is worth noting also 
that the two modifications proposed in sections [2] and El could be combined to obtain a 
stable 2-TL scheme together with thermal NL residuals of smaller magnitude. 

More generally, the aim of this paper is to emphasize that there may be a considerable 
benefit to relax the unnecessarily constraining principles (i)-(ii) for the design of all kinds 
of implicit schemes based on a linear separation (SI and ICI schemes). In practice, starting 
from an initially unstable scheme obtained through the traditional approach, a dramatic 
improvement may sometimes be obtained if the approach proposed here is used for modi- 
fying this initial scheme in only slight details. The discussions in this paper clearly do not 
offer a complete picture of the properties of the modified schemes proposed above com- 
pared to their traditional counterparts. Before extending such modifications to the actual 
NWP framework, it would be necessary to evaluate more precisely their practical impact 
in terms of accuracy and response corectness on forecasts performed in real conditions. 
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List of Figures 

Fig. 1 : Growth-rate T with the HPE system as a function of the the tropopause 
location a for the thermal profiles and discretisation settings indicated in the text. Circles : 
traditional 3-TL SI scheme ; crosses : modified 3-TL SI scheme. 

Fig. 2 : Growth-rate T as a function of «, in the limit of long time-steps for the complete 
EE system examined in section O Solid line : traditional 2-TL SI scheme ; dotted line : 
modified 2-TL SI scheme. 

Fig. 3 : Evolution of the spectral norm of vorticity £ and horizontal divergence D for 
a real-case with a 2-TL SI EE system. Solid line : vorticity (right axis) ; dashed line : 
divergence (left axis) ; thin line : traditional 2-TL SI scheme ; thick line : modified 2-TL 
SI scheme. 
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Fig. 1 - Growth-rate T with the HPE system as a function of the the tropopause location 
a for the thermal profiles and discretisation settings as indicated in the text. Circles : 
traditional 3-TL SI scheme ; crosses : modified 3-TL SI scheme. 
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Fig. 2 - Analytical growth-rate T for the EE system as a function of a in the limit of 
large time-steps. Solid line : traditional 2-TL SI scheme ; dotted line : modified 2-TL SI 
scheme. 



20 



1.5 



0.5 




1.5 



0.5 



t(h) 



FlG. 3 - Evolution of the spectral norm of vorticity ( and horizontal divergence D for 
a real-case with a 2-TL SI EE system. Solid line : vorticity (right axis) ; dashed line : 
divergence (left axis) ; thin line : traditional 2-TL SI scheme ; thick line : modified 2-TL 
SI scheme. 
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